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Abstract 


We investigate the inner crust structure of neutron stars using the Skyrme- 
Hartree-Fock approach with the Coulomb interaction treated beyond the 
Wigner-Seitz approximation. Our results suggest that the shell effects as¬ 
sociated with unbound neutrons play an important role and, in particular, 
lead to complicated phase transition pattern between various nuclear phases 
(as a function of the density). Namely, we show that the relative energies of 
different phases are rapidly oscillating functions of the neutron density. In 
the semiclassical approach this behavior is explained as an interference effect 
due to periodic orbits of similar lengths. We discuss also the dependence of 
the shell effects on pairing correlations. 
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I. INTRODUCTION 


The outer layer of a neutron star has a complex structure which depends strongly on the 
nuclear density. In the inner crust of the star, due to the high density and pressure, a large 
fraction of neutrons occupies unbound states. Nuclei which are still present are therefore 
immersed in a neutron gas. The structure of this part of the star has been the subject of 
a considerable theoretical effort MU- The main conclusion emerging from these studies is 
that the nuclei form exotic phases, stabilized by the Coulomb interaction, and characterized 
by various crystal lattice structures. 

The theoretical predictions on the structure and stability of the nuclear phases are based 
on the minimization of the energy of the neutron-proton-electron (npe) matter within a 
Wigner-Seitz cell, whose radius is treated as a variational parameter. An agreement has been 
reached on the existence of five phases formed in the region where the nucleon density varies 
from 0.03 to 0.1 fm -3 . The following chain of phase transitions is predicted as the density 
increases: spherical nuclei —> rods —> slabs —> tubes —> bubbles —> uniform matter. The 
physical mechanism responsible for the phase transition pattern results from the interplay 
between the Coulomb and surface energies. However, the differences between the energy 
density of the various phases are only of the order of a few keV/ fm 3 , for a wide range of 
densities |T3[. The precise transition pattern is therefore dependent on the ingredients of the 


theoretical models. For example, recent calculations based on the self-consistent Hartree- 
Fock (HF) method have predicted that the phase composed of spherical nuclei is stable up 
to a density around 0.077 fm -3 , above which the uniform matter is energetically favored 
17H ■ Other phases appear at higher energies and form excited configurations. 

In the present paper we focus on an effect which has been overlooked in previous studies 
but should strongly influence the phase transition pattern. In all to date approaches, shell 
effects have been either completely neglected, as in models based on a liquid drop formula 
or in Thomas-Fermi models, or have been taken into account only for bound nucleons. On 
the other hand, it has been recently shown in the semiclassical approximation, based on the 
Gutzwiller trace formula, that the shell effects in inhomogeneous nuclear matter are quite 
strong and yield corrections to the energy density of the order of several keV/fm 3 . This 
value is comparable to the differences between the liquid drop energy densities predicted for 
the different nuclear phases ||TB|-[?(]|j . Moreover, the shell energy oscillates as a function of 
density and depends strongly on the spatial order in the system. In particular, its amplitude 
depends both on the nuclear geometry and on the lattice type. 

The shell energy of npe matter in the neutron star crust has two origins. The first 
one is associated with bound nucleons and has been shown to be small by Oyamatsu and 
Yamada in Ref. |n|. The second one is related to unbound neutrons which may scatter 
on inhomogeneities and form resonant states. The situation is somewhat similar as for the 
Casimir effect in quantum field theory and condensed matter (see e.g. Refs. ]2T| fT]|| and 
references therein), where, in some systems, a fluctuation induced interaction leads to a 
correction to the total energy. In the same way, the scattering of unbound neutrons lead to 
an effective interaction between nuclei immersed in the neutron environment. This effect can 
be termed the fermionic Casimir effect and has been studied for simple geometries in Refs. 
mmmm- For impurities that are small compared to the neutron Fermi wavelength, the 
generated interaction is mainly due to s-wave scattering. It gives then rise to the Ruderman- 
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Kittel correction which is small compared to the Fermi energy of neutrons |f2d|.p?7 


However, the s-wave scattering limit is not valid in the neutron star crust, since the 
neutron Fermi momentum kp is of the order of 1/m -1 and thus kpR > 1, where R is the size 
of an inhomogeneity. Consequently the contribution of higher partial waves to the scattering 
process cannot be neglected, and gives rise to much larger energy correction, as for large 
objects more of the incident wave will be reflected ||T8|- ^0|^ ,^5|] . The leading shell energy 
contribution, associated with the neutron scattering, to the total energy can be determined 
in the semiclassical approximation for various geometries of nuclear shapes (see appendix B). 
It has been shown to yield energy corrections that influence the structure of inhomogeneous 
nuclear matter | [T8| |20|| . It is also obvious that since the gross shell structure is determined 
by the shortest periodic orbits in the system, the shell energy is a sensitive function of the 
crystal lattice type, as well as of the shape of inhomogeneities [|18|,^5j. Unfortunately the 


semiclassical approach has drawbacks. First, it assumes that the obstacles are impenetrable 
scatterers which may overestimate the amplitude of the shell energy. Second, the method 
lacks of the mutual interplay between the shell energy and the liquid drop energy. Such a 
coupling is quite important since it allows a part of the shell energy to be “absorbed” into 
deformation. Hence a microscopic treatment of the problem is needed where both effects, 
i.e. the one coming from the liquid drop energy and the shell energy term are treated on 
the same footing. 


II. THE METHOD 


The ground and excited (isomeric) states of npe matter have been determined by the 
Skyrme-Hartree-Fock method. We have solved the HF equation by discretization in coor¬ 
dinate space, within a rectangular box and with periodic boundary conditions imposed on 
the nucleon wave functions |28| . This technique has the advantage that it allows to describe 


any kind of shapes of the nuclear density distribution with the same accuracy, and to in¬ 
troduce naturally periodic boundary conditions. The code that was previously developed to 
study nuclei with an ellipsoidal symmetry has been modified to accommodate the periodic 
solutions. 

The high efficiency of the method can be attributed to the iterative procedure used to 
solve the HF equations, the imaginary time step method p9[. One of its features is that only 


those single-particle wave functions which contribute to the mean field are calculated. It is 
particularly important in our case, since the densities which we are dealing with here require 
to consider up to 1000 nucleons per box. In order to take them into account properly, we 
have used 600 wave functions. Due to Kramers degeneracy, they correspond to 1200 single¬ 
particle states and thus form a sufficiently large Hilbert space for our purposes. 

As a starting point of the HF iterations, Fermi gas wave functions have been used, which 
turned out to provide a good first guess of the true ground state. The criterion for the 
termination of the iteration process is the total energy difference obtained after a fixed 
number of imaginary time steps. In our calculations we terminated the iteration procedure 
when the differences between the total energy calculated at two successive iterations did not 
exceed 1 keV. For stable configurations this value was of an order of magnitude lower. 

An important parameter of a mesh is the spacing between the mesh points. It must 
be small enough to lead to accurate results and large enough to accommodate the very 
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large number of wave functions needed to describe the neutron gas. We have chosen as a 
compromise between these two requirements a value of the mesh size equal to Ax = 1.3 fm. 

This choice of Ax fixes a bound to the maximum momentum of a nucleon that can 
be represented by the discretized wave function. The highest momentum is indeed equal 
to k cr it = 7t/Ax ~ 2.42/m , and corresponds to the maximum kinetic energy e cr it = 


121 MeV. 


2mAx 2 

There are two main sources of errors caused by the discretization. One of them is related 
to unbound neutrons which may occupy high energy states and thus can not be properly 
described by a discretized wave function. However even at the highest neutron density that 
we are interested in here, i.e. p n ~ 0.1 fm -3 , the neutron Fermi momentum is equal to 
kp 1.44 fm -1 which is far below the critical value. A second source of errors is associated 
with nucleons that are bound in a nucleus. Typically, a nuclear potential at densities p > 0.03 
fm -3 has a depth which does not exceed 20 MeV and a radius lower than 10 fm 


12 . In a 


Fermi gas approximation for bound nucleons, the momentum of the last bound nucleon is 
around 1 fm -1 , well below the maximal value that can be represented on the mesh. Let us 
note that the Fermi gas approximation gives an upper bound of the maximum momentum 
since the wave function of a bound state has a long tail due to the large diffuseness of a 
nuclear potential immersed in a neutron gas. 

To decrease the numerical task, we have assumed that the nuclear density distribution is 
symmetrical with respect to three planes, allowing however the existence of triaxial shapes. 
Parity breaking configurations are thus excluded. 

Different phases may coexist for a given nuclear density. To determine them and to study 
their stability requires the construction of adiabatic paths between phases. This has been 
done by introducing a constraint on the quadrupole moment of the proton density distribu¬ 
tion in the Skyrme HF equations. In order to compare energies of configurations differing 
by the proton quadrupole moment, a careful treatment of the Coulomb interaction, both 
between protons in different nuclei of the lattice, as well as between protons and electrons, 
is needed. It has lead us to solve the Poisson equation for the electric potential within a 
box with periodic boundary conditions. This approach is clearly free from limitations of the 
usual Wigner-Seitz approximation, where the influence of the lattice type on the Coulomb 
energy is not taken into account. The electrons are assumed to form a uniform relativistic 
gas where the screening effects are neglected. The total charge of the box is zero [IT| , which 
gives a relation between the electron and proton densities. 

Summarizing, we solve the Hartree-Fock equations for the npe matter with contributions 
in the energy functional coming from the nuclei, the neutron gas and the electrons. In such 
a calculation, the liquid drop and the shell energy parts of the energy are automatically and 
self-consistently included. The total energy of the system can be expressed in the form: 

EtotfyPpi Pm Pe ) A 'Skyrme^yPpi Pn) T T'Comz(Pp) Pe) T E e i ( p e ), (1) 

where E S k y rme comes from the Skyrme energy functional |2l| , E Cou i is the Coulomb energy 
term, and E ei denotes the electron kinetic energy term. The densities p n , p p , p e denote the 
average densities of neutrons, protons and electrons, respectively. They are related to the 
density distribution through the equation: 

Pi = y J v Pi(r)d 3 r, i=n,p,e, (2) 
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where V is the volume of the cell. The charge neutrality condition implies that p e = p p . 
The total nuclear density is given by p = p n + p p . In the calculations we have used the SLy4 


Skyrrne force. This recent parametrization [|30| gives a good description of neutron matter 
at low densities pT|,|T7|. For a few test cases, calculations have also been performed with the 


SLy7 force which does not lead to any qualitative differences. 

The shell energy which is implicitly included in the expression (HD can be decomposed 
into two parts. The first one, denoted by E '^ eH , is determined by the density of neutrons 
p™ and protons p* n inside the nucleus i.e.: 




VniLCl JVnucl 


Pi(r)d 3 r , 


n ,P, 


( 3 ) 


where V nuc i is the volume of the nucleus, 
neutrons p° ut and protons p° ut , where 


The second, E^. u , is a function of unbound 


P° Ut = 


V - V nud JV-V n . 


Pi(r)d 3 r, 


!=n,p. 


( 4 ) 
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17], the contribution to the energy due to 


Since p^ is practically zero for p < 0.08 fm 
unbound protons is marginal. 

It should be emphasized that the decomposition of the density into an inner and an outer 
part is somewhat arbitrary. There are various methods of defining these quantities that have 
been used in liquid drop based approaches |L5],[T6|. In the HF method this decomposition 
is not necessary to calculate the total energy of the system. Nevertheless, in order to 
understand its behavior as a function of physical parameters such as nuclear density, proton 
fraction, lattice constant, etc., it is important to identify and analyze the origin of the various 
contributions to the total energy. 


III. RESULTS 


The HF equations are solved for a fixed total density p, box size d and with a constraint 
imposed on the proton quadrupole moment. In this way, it is possible to determine not 
only the ground state but also excited configurations, differing from the ground state by, for 
example, the lattice type. Since we solve the problem in a cubic box with three symmetry 
planes, several lattice geometries can be generated like e.g. simple cubic crystal (see), face 
centered crystal (fee), or body centered crystal (bcc) p2| . 

In order to verify that the absolute minimum of the functional ([T) has been found, we 
have performed a set of constrained calculations corresponding to various values of the proton 
quadrupole moment for each value of the total nuclear density. It should be emphasized that 
we did not perform the total energy minimization with respect to the parameters of the box 
which were kept fixed at values corresponding to those obtained in the reference jl2| . 

In the figure 1 we show the ground state energy density as a function of the proton 
fraction: Z/A = p p /p (Z = Vp p and A = V p, where V = d 3 is the volume of the box) for 
four values of the nucleon density. The calculations indicate that for the density range we 
consider here, the optimal Z/A ratio varies between 0.03 and 0.04. 

For a total density p = 0.055 fm -3 , the energy density is plotted in Fig. 2 (middle 
part of the left panel) as a function of the proton quadrupole moment Q p = Q 2 0 . Most 
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studies predict that, at this density, the nuclear spherical phase is energetically favored (see 
e.g. [|l|,|T^] and references therein). One can see that in our self-consistent calculations, the 
energy minimum appears for a non-zero quadrupole moment. In order to identify the shape 
of the nuclei, as well as the lattice type, we have plotted the neutron and proton density 
distributions integrated along the z and y axes: 


f d 

Pi{%,y) = / pi{x,y,z)dz, 

(5) 

Jo 

rd 


Pi( x , z )= Pi( x ,y, z )dy, 

(6) 

Jo 

(7) 


where i = n,p. In the lower and upper parts of the left panel of Fig. 2, the integrated 
density distributions are plotted for both protons and neutrons. We have shown the density 
distributions for the configurations indicated by arrows, which correspond to the spherical 
and the most elongated nuclear shapes. The density distributions in the energy minimum 
(not plotted) have a pattern very similar to the spherical one. They correspond to a see 
lattice of deformed nuclei immersed in a neutron gas. 

An excited minimum appears for a proton quadrupole moment of about 2500/m 2 . It is 
associated with a rod-like (“spaghetti”) phase (upper part of the left panel of Fig. 2). The 
barrier between these minima is of the order of 0.5 keV/ fm 3 . One can see that contrary to 
previous calculations, the spherical phase is no longer the most favored one and appears at 
an excitation energy of about 1 keV/fm 3 , even above the “spaghetti” phase. When keeping 
the proton number fixed and adding neutrons to the system, the energy density pattern 
slightly changes. The right panel of Fig. 2 corresponds to the situation when the number 
of neutrons have been increased to reach the total density p = 0.057 fm -3 . As shown in 
the right panel of Fig. 2, the position of the minimum is not altered. It suggests that the 
deformation of the ground state is mainly related to the bound nucleons. On the other 
hand, there is a substantial change of the position of the spherical phase which drops down 
in energy and is now below the “spaghetti” phase. 

The Fig. 2 shows that the change in the neutron density does not influence qualitatively 
the structure of the ground state since the energy minimum stays approximately at the same 
proton quadrupole deformation. However it is not always the case as one can see from Fig. 
3, where the energy density has been plotted as a function of proton quadrupole moment 
for larger total densities. The left panel shows the results for p = 0.075 fm -3 , and the right 
one for p = 0.079 fm -3 . The two almost degenerate minima visible in the left panel differ 
by the lattice type only, and both correspond to nuclei immersed in a neutron gas. The 
spherical configuration corresponds to a see lattice and the deformed one to a bcc lattice. 
There is also an excited configuration associated with a slab-like (“lasagna”) phase. Still 
increasing the number of neutrons, but at fixed proton number, one sees on the right panel 
of Fig. 3, that the relative position of the minima is not stable with respect to the neutron 
density. The spherical phase is no longer energetically favored and the bcc configuration 
merges with the slab phase, forming the lowest energy configuration in the form of “rippled” 
nuclear slabs. 

The question arises whether these changes are related to the shell effects or to the liquid 
drop part of the energy. In order to better visualize this effect, we have plotted (see Fig.4) 
the energy differences between different configurations: 
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• the spherical nuclear phase (Q p = 0 fm 2 ), 


• the “spaghetti” phase ( Q p = 2400 frn 2 ), 

• the “lasagna” phase ( Q p = 2600 frn 2 ). 

The solid curve represents the quantity A E/V = ( E spherica i — E spaghetti )/V while the dotted 
curve corresponds to A E/V = ( E sp h er icai — Ei asagna )/V . One notices that these differences 
oscillate rapidly as a function of the total density. It seems therefore unlikely that the liquid 
drop part of the energy is responsible for such rapid fluctuations. On the other hand, the 
leading shell energy contribution to the total energy, for two obstacles located at the distance 
d and immersed in the neutron gas, can be determined in the semiclassical approximation 
(see appendix B): 
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where m n is the neutron mass and i — 0,1,2 stands for two spherical, cylindrical and planar 
obstacles, respectively (we assume that rods and slabs are parallel to each other). In the 
above equation, L defines the length of the obstacle, and R is its radius (in the case of a 
slab it is defined as half of its width). This expression exhibits oscillations, where the gross 
staggering pattern comes from the shortest periodic orbit of length d between nearest nuclei. 
This can be used to estimate qualitatively how the difference between the shell energies of 
the two minima varies as a function of the neutron density. 

Assuming that the distance between the nearest nuclei is given by r, the shell energy 
contribution clue to unbound neutrons has the form: 


£Soccos(2*;M. (9) 

The quantity plotted in Fig. 4 represents the difference between the energy density of two 
phases, and if it is associated with shell effects, its behavior should be expressible as a sum 
of cosine (or sine) functions. In particular, the leading contribution with the longest period 
comes from the single repetition of the shortest periodic orbits: 

^ E °heii « A COS (2fc£Vi) + Aj cos(2fc£Vj) i ± j, (10) 

where i,j = 0,1,2 for spherical nuclear phase, the “spaghetti” phase, and the “lasagna” 
phase, respectively. The coefficient Ai is related to the stability of the orbit, , r l is the 
distance between nuclear inhomogeneities, and k™ stands for the Fermi momentum in the 
phase i. One can infer from this expression that the interference between the contributions 
of different orbits is at the origin of the fluctuations of AE° p * u as a function of the density. 
Indeed, taking k™ — k v F ~ (37r 2 p" u *) 1 / 3 and assuming that A { « Aj = A we get: 

AE °sheii ~ 2 A sin {2k n F rij) sin(^Ar„), (11) 

where = 1 (r t + r 3 ) and A= (r t — r 3 ) . Hence the shell energy difference oscillates as a 

7 r 

function of the Fermi momentum with a period Akp = —, or equivalently, as a function of 

Aj 

h 2 n 2 

the neutron Fermi energy, with the period A €p = - 

2 m n rtj 
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This quantity can be estimated from our HF calculations. Let us define the radius of the 
nuclear inhomogeneity as the distance between the center of a nucleus and the point where 
the neutron density decreases to the value ^(p° ut + p™). This gives a radius R 0 ~ 6.23 fm 
for the spherical and R\ ~ 3.86 fm for the “spaghetti’' phases, at a density p ~ 0.055 fm _3 . 
Since the length of the lattice d is in this case equal to 26/m, the distance between the two 
nearest spherical-like and spaghetti-like objects can be estimated as tq = d—2R 0 « 13.54/m, 
and r\ = d — 2Ri « 18.28/m, respectively. This gives a period of oscillation equal to 
Ae/ ~ 0.8 MeV. On the other hand, the HF calculations predict the neutron Fermi energies 
to be equal to: 


e n F {Q p = 0 ,p = 0.063) = 12.043MeV, 
e n F (Q p = 0 ,p = 0.053) = 10.74 MeV, 
e n F (Q p = 2400, p = 0.063) = 11.931 MeV, (12) 

ef(Q p = 2400, p = 0.053) = 10.607MeV, 

which give the periods: 

A e n F {Q p = 0) = 1.303 MeV, (13) 

Ae F (Q p = 2400) = 1.324 MeV. (14) 

They have to be compared with the 0.8 MeV obtained on the basis of the semiclassical 
analysis. The difference between these values is mainly related to the fact that in the 
semiclassical approximation we have neglected both the coupling of the shell effects to the 
liquid drop terms, and higher order corrections related to other orbits (including the multiple 
repetitions of the shortest orbits). Nevertheless, taking into account how drastic these 
approximations are, the agreement is surprisingly good. It is even better for the case of 
the energy density difference between the spherical and the slab phases (dotted curve in 
the Fig. 4). The radii are then: R 0 « 5.9/m, and i? 2 ~ 3.59 fm at a density p ~ 0.075 
/m -3 , with the “radius” of the slab defined as half of its width. The distances between the 
inhomogeneities for these two phases are ro ~ 9/m and r 2 ~ 13.62/m, respectively. Half of 
the oscillation period is equal to \Aef ~ 0.801MeV. On the other hand, according to the 
HF calculations the neutron Fermi energies are equal to: 

e F (Q p = 0 ,p = 0.079) = 12.645MeV, 
e r f(Q p = 0 ,p= 0.0749) = 11.779MeV, 
e n p(Q p = 2600, p = 0.079) = 12.759 MeV, (15) 

e F (Q p = 2600, p = 0.0749) = 11.884MeV, 

which give for half of the periods the values: 

^A 6 n F (Q p = 0) = 0.866MeV, (16) 

\Ae n F (Q v = 2400) = 0.875 MeV. (17) 

(18) 
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The agreement with the semiclassical estimation is very good. This demonstrates that the 
oscillatory effect of the energy density difference between the phases is due to the shell 
energy of unbound neutrons. From eq. (|TT|), it is clear that another oscillatory effect may 


appear with a much longer period Akp = 


2n 
A n 


Its mechanism is somewhat similar to the 


well known supershcll phenomenon in atomic clusters [Q. This effect however, requires the 
stability of one phase over a wide density range and therefore may be not well pronounced. 

Up to now, we have neglected the pairing interaction. The pairing energy is a smooth 
function of the density and is not expected to play an important role for the liquid drop 
part of the total energy expansion. On the contrary, pairing correlations may considerably 


influence the shell energy fl34|j35 


The pairing energy can be divided with a good approximation into two parts. The 
first one comes from the pairing interaction acting among bound nucleons, and the second 
one is associated with pairing of unbound neutrons. These two contributions can be, in 
the first approximation, treated as independent quantities, because of the short range of 
the pairing interaction. For example, in the usual zero-range approximation of the pairing 
interaction (delta force) the amplitude for pair scattering is proportional to the integral: 

J 0^(r)0 z *(r)0 m (r)0 n (r)(i 3 r, where </>&(r) are the single-particle wave functions. On the other 

hand, bound single-particle states, as well as resonant states inside nuclei have a small overlap 
with states forming continuum. Hence the integral describing the scattering of a Cooper 
pair from nuclei to continuum or vice versa will be relatively small, and can be, in a first 
approximation, neglected. 

The pairing interaction acting among bound nucleons will decrease the amplitude of 
shell effects inside the nuclei. Note however that these effects can at most lead to a slight 
deformation change of nuclei and do not influence the phase transition pattern. Hence in 
the following we focus only on the pairing energy associated with the neutron gas. Since the 
pairing gap A for neutron matter at subnuclear densities reaches 3 — 4 MeV at the Fermi 
level |3B|] , the occupation of the single-particle states close to the Fermi level are significantly 
modified by pairing. In order to analyze the influence of pairing on the shell structure, let 
us consider a Fermi gas of neutrons in a volume Fat a density interacting via a pairing 
force. The total energy of the system can be expressed in the BCS approximation as: 


E = 


/ OO 

^l(e)g(e)de + E 1 

-OO 


pairi 


(19) 


where g(e) is the density of single-particle states at the energy e and /i is the chemical 
potential defined by the condition: 

1 


r oo 

Pn t =yj_ 


The coefficient v 2 (e) represents the BCS occupation factor: 


(0 = 


1 - 




(e - g,) 2 + A(e) 2 y 


( 20 ) 


( 21 ) 


where A(e) is the (energy dependent) pairing gap. In order to avoid divergences, we assume 
that the function A approaches zero sufficiently rapidly to make the integral in equation 
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(|19|) convergent. The term E pair denotes the pairing energy. If the Fermi gas is subject 
to boundary conditions, e.g. due to presence of impurities in the system, the level density 
contains a fluctuating term: 


d) =j(e)+9c(«,d), 


( 22 ) 


where g(e) is the smooth part of the level density. The fluctuating part gc (e, d) depends on 
the shape and the mutual geometrical arrangement of the impurities parametrized by d and 
is responsible for the shell correction energy. For two spherical, cylindrical or planar obstacles 
located at a distance d, the semiclassical approximation gives the leading contribution to gc 
in the form (see appendix B): 


1 


f/i) ‘ e'"cos(2 kd-ij), 


(23) 


where e = 


h 2 k 2 
2 m n 


and i takes the values 0,1,2 for the various shapes (compare to eq. 


The total energy of the system can now be divided into smooth and oscillating parts (see 
appendix A): 


E = E + E c + A E + E, 

roo 

= / wl(e)g(e)de + 


pair i 
oo 
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r oo 
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ev?(e)gc(e, d)de 


v U e )g c ( e ’ d ) de - 


5-oo^l{e)~g{e)de 




+ E pa i r + 0 ( 1 / 1 /), 


(24) 


where the first term represents the smooth part of the total energy. The term Ec is the 
Casimir energy, and A E provides a correction to the Casimir energy due to the Fermi 
energy renormalization in the presence of obstacles ||TH| , P5|| . The smooth part of the chemical 
potential fi is calculated from the condition: 


N = 




(25) 


It implies that the neutron shell energy in the presence of pairing can be expressed in the 
form: 


E2U = E c + AE + 0(1/V% 


(26) 


where we have neglected corrections due to the change of pairing energy in the presence of 
obstacles (see appendix A). 

In Fig. 5 we show the behavior of the shell energy as a function of the distance between 
obstacles without pairing and with a pairing interaction characterized by a pairing gap 
equal to A = 3.5 MeV for k < 20/m _1 (a cutoff has been introduced to avoid divergences). 
This pairing strength corresponds to the pairing gap at the Fermi level for a neutron Fermi 
momentum kp « 1 fm _1 obtained with the Argonne V\± and Uis interactions [[36). One can 


see that the amplitude of the shell energies are smaller but the decrease is lower than 10%, 
as compared to the no pairing regime. We conclude from this comparison that the pairing 
interaction will not change qualitatively the HF results. 
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IV. CONCLUSIONS 


In this paper we have analyzed the structure of the neutron star crust at zero temperature 
using the Hartree-Fock approach with an effective Skyrme interaction. In particular, we 
focused on the problem of shell effects associated with unbound neutrons scattered by nuclear 
inhomogeneities. Our results show that: 

• the relative energies of different phases fluctuate rapidly as a function of the total 
density, 

• the fluctuations can be attributed to the shell effects associated with unbound neu¬ 
trons, 

• the neutron shell energy density is of the same order as the liquid drop energy density 
differences between the various phases present in the crust, and its behavior can be 
easily understood in terms of the periodic orbits in the system, 

• the pairing correlations lead only to a slight decrease of the shell energy and are not 
expected to influence significantly the structure of the crust. 


Consequently, we may conclude that the structure of the crust may be quite complicated 
since the shell effects associated with unbound neutrons may easily reverse the phase tran¬ 
sition order predicted by liquid drop based approaches. Moreover the number of phase 
transitions may increase since the same phase may appear for various density ranges (see 
Fig. 4). 

Our results suggest also that several phases, differing by the lattice type, could coexist 
in the crust. This possibility was not taken into account in previous investigations since in 
a liquid drop based approach, the system favors only one lattice type for a given nuclear 
shape. However, since the neutron shell energy is very sensitive to the spatial order in 
the system, this simple picture will change, and one may expect the coexistence of various 
lattice geometries with different lattice constants. Note also that once a phase is formed, 
its geometric order will be stabilized by both the Coulomb repulsion and the shell energy. 
Indeed, although the Coulomb energy is a smooth function of the nuclear displacement, the 
shell energy is not. Since several different orbits contribute to the shell effects (except for 
the slab-like phase) the displacement of a single nucleus from its equilibrium position in the 
lattice will give rise to interference effects depending on the lattice type. Hence it is likely 
that the system will favor distorted lattices, or lattices with defects which decrease the shell 
energy. In that respect the present results support the conclusions drawn on the basis of a 
simple model 


All these arguments suggest that the inner crust is an extremely rich and complicated 
system which may turn out to be on the verge of a disordered phase. 
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APPENDIX A: SHELL ENERGY OF THE FERMI GAS WITH IMPURITIES 


We derive here expressions for the shell energy associated with impurities immersed in 
a fermionic environment. 

Let us consider a Fermi gas confined in a volume V, large compared to both the Fermi 
wavelength (kpV » 1) and the volume of the obstacles. At the end of the calculations the 
limit V —> oo will be taken. The shell energy is an interaction energy between “foreign” 
objects immersed in the gas. Its origin comes from the modification of the level density by 
the presence of obstacles. The total level density can be decomposed in the following way 




0(e,d) = g(e) + g c (e, d), 


(Al) 


where g is the smooth part of g and can be obtained either from the Weyl expansion [37.33 
(in the case of impenetrable objects) or within the Thomas-Fermi approach [[33]]. The set of 
parameters d = (<A, d,2, •••, d n ) describes the geometrical arrangements of the obstacles. The 
smooth part can also be treated as the total density of states for a system with obstacles 
located at infinity. Namely: 


lim g(e, d) = g(e), for i=l,2,...n. (A2) 

CL'i ^ oo 


The fluctuating part of the level density gc is associated with the formation of resonances. 
Its evaluation in a semiclassical approximation is presented in the next section. 


The shell energy can be defined in a form similar to the Casimir energy [IS-20,24,25,32 


ri‘ rn 

E ahe ii( d) = / deeg(e, d) — / deeg(e), 

J— oo J— OO 


(A3) 


where g and g denote the chemical potentials associated with the total and smooth level den¬ 
sities, respectively. They are related to each other through the particle number conservation 
condition: 


N 


rM |'/J 

/ deg(e, d) = / deg(e). 

J —oo J—oo 


(A4) 


The decomposition of the integrals in equation (|A3| ) results in the following expression: 

E sh eiiiA) = AE + E C = [ eg(e)de + f eg c (e,d)de, (A5) 

J LL J —OO 


where A E represents the correction to the shell energy due to the variation of the chemical 
potential induced by the presence of obstacles, while Ec corresponds to the ’’pure Casimir- 
like” force. Since Ag = g — g oc 0(1/V) and g oc V then 
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eg{e)de = gA gg(g) + 0( 1/V), 


g(e)de = Agg(g) + 0{ 1/V), 


(A6) 


and thus it is clear that both A E and Ec do not vanish in the limit of infinite volume. The 
condition (|A4|) can be used to calculate A g: 


Ag = 


<?0) 


gc(e,d)de + 0(l/V 2 


The shell energy can be finally written in the form: 


J shell 


(d) = / (e-g)g c (e,d)de + 0(l/V). 

J — OO 


(AT) 


(AS) 


Taking the limit V —» oo, the higher order terms vanish and the hrst term in the above 
formula yields an exact expression for the shell energy. 

In the case of a superfluid Fermi gas, the energy gap depends on the momentum k. 
To avoid divergences in the expression for the total energy, one has to assume that A (k) is 
tending to zero at least like 1 /k for k —> oo. The pairing properties of the system are affected 
by the presence of obstacles and thus the solution of the gap equation with impurities A (k, d) 
differs slightly from the smooth pairing gap A (k) obtained using the smooth level density: 
A (k, d) = A (k) + A c (k, d). A c denotes the quantum correction to the smooth pairing gap 
and is due to the fact that the level density is modified. Hence, eq. (|A3|) is generalized for 
superfluid systems in the form: 


E she ii{ d) = / deev 2 (e,d)g(e,d) 


deevl(e)g(e) 


Epair (d) - E 


pair i 


(A9) 


where v 2 is the BCS occupation factor in the presence of obstacles (|2l|) and the term vp, 
represents the smooth occupation factor obtained for impurities far apart from each other. 
The last two terms are due to the correction of the pairing energy related to the appearance 
of resonant states between obstacles. 

For simplicity, we assume that A c = 0 (and thus v 2 = v 2 ) and neglect the effect of 
the geometrical arrangements of the obstacles on the pairing properties of the system. The 
equation simplifies then to the form: 


Eshell (d) — 


deevl(e)g(e), 


deev^g^d) - i 

J—oo J—oo 

where the chemical potentials g, and g are determined from the conditions: 

/ oo r oo 

devl(c)g{e, d) = / dev 2 p(e)g{e). 

-OO J—OO 

Expanding v 2 around g one gets: 

dv 2 p(e) 


(A10) 


(All) 


Esheiii.d') A g 


f 00 dee-E^- 

l—oo dfi 


/ OO / ( 

^ dee ( v 2 p(e) + Ag- 


dg 


gc(e,d) + 0(l/V). 

(A12) 
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Note however that the term: Afi 




dfi 


gc(e,d) is also proportional to 1/V and 


thus vanishes when V —> oo. Determining A fi from condition (|A11|) one gets the expression 
for the shell energy: 


J shell \ 


(d) = / evl(e)g c {e,d)de- vl(e)g c (e,d)de 


poo Q 2 

f-oo e ~dp, v p. 


vl(e)g(e)de 


J-oc f,vl(e)~g(e)de 


+ 0(1/D). (A13) 


Since v~(e) —> 9(g — e) and 


dvUe) 


' dp, 

( |A13|) becomes identical with (|A3| ) in the no pairing limit. 


— 5(e — fi) when A (k) tends to zero, the equation 


APPENDIX B: EVALUATION OF THE FLUCTUATING PART OF LEVEL 
DENSITY AND THE SHELL ENERGY IN THE SEMICLASSICAL 

APPROXIMATION 


We derive here expressions for the shell (interaction) energy for two identical, impene¬ 
trable obstacles at large separations (||). Three geometries of the obstacles are considered: 
spherical, rod-like and planar geometries. The system is characterized by a set of parameters 
i?, L, d which denote the radius, length (in the case of rods and slabs), and the distance be¬ 
tween obstacles, respectively. The distance is measured as half of the length of the shortest 
periodic orbit. It means that e.g. in the case of spheres the distance between their centers 
is equal to d + 2 R. In the case of a slab, R denotes half of its width. 

In the semiclassical approximation, the fluctuating part of the level density gc can be 
obtained through the summation of the contributions associated with the classically allowed 
periodic orbits in the system. Recent exact calculations based on the S-matrix approach 
H indicated that this approximation work surprisingly well even if the obstacles are close 


to each other, i.e. when d « R. In general, the summation is a tedious task if many orbits 
contribute to the level density. Especially when elliptic or parabolic orbits are considered, 
the problem requires a more careful treatment in order to avoid divergences j33| . The 


semiclassical approximation to the interaction energy between two obstacles is particularly 
useful if the periodic orbits of the system are unstable (hyperbolic) and isolated [^3[. It 
holds true for spherical obstacles. For two spheres of radius R, located at a distance d, the 
contribution to the level density arising from the single periodic orbit is: 


d) = 


md 


E 


Tifik ^ 2 skill (2/ 


7T 

^cos(2 nkd-a-), 


(Bl) 


where A = 2 log 




is the Lyapunov exponent determining the sta¬ 


bility of the orbit, cr is the Maslov index depending on the type of boundary conditions, and 
h 2 k 2 

e = ——. For Dirichlet boundary conditions, which we assume here, each reflection changes 

the phase by 7r, thus in our case a = An where n is the number of repetitions of periodic 
orbit. The shell energy calculated from the expression (|A3|) gives ||18||: 
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rp sphere 
^ shell 


(d) 


h 2 

167T md 2 


OO 


E 

n =1 


1 

n 3 skill 2 (ky) 


(2 nk F d cos(2 nk F d) 


sin(2 nkpd)). 


(B2) 


At large distances, i.e. for — >> 1, the sum can be truncated at n — 1 and the expression 

R 

simplihes to: 


jp sphere 
^ shell 


(d) 


n 2 R 2 
8 m 



1/3 1/3 

ff 


cos(2 k F d) + 0(l/d 5 ), 


(S3) 


where we have expressed the amplitude of the shell energy in terms of the density p 
This expression is associated with the fluctuating level density: 


p 

Ay J? 

37r 2 


d) « 


1 mR 2 

2 7T h 2 kd 


cos(2 kd) + 0(l/d 3 ). 


(B4) 


In the case of two parallel rods of length L, the orbit is no longer isolated. Hence it is 
convenient to consider first the two-dimensional case of two disks. Then the single periodic 
orbit gives rise to the fluctuating part of the level density: 


9c s \k,d) 


md ^ 1 

7T h 2 k n=1 sinh( 2 y) 


cos(2 nkd). 


(B5) 


In order to obtain the g™ d one has to calculate first the fluctuating part of the particle 
number iV™ d (/i, d) = [ g r c d ( e, d)de. But N™ d has to be an integral of the product of 

gfj sk an( ] ^-j ie level density of one-dimensional Fermi gas. Namely, the quantity Nq oc1 can be 
expressed as: 


N r c ° d (k F ,d ) = 


27t Jo 
L ki 


k f f Ti 2 k 2, 

g d s sk m 




E 


47T sinh(^) 


2m J J-^/k 2 F -k 2 

Ji(2nk F d), 


dkr. 


(B6) 


where k p denotes the momenta parallel to rods and J\ is the Bessel function of the first kind. 


v p 

fl /V 777 

Since g r c od {k,d ) = — we getfl: 


dN r c od 


9 r c od fd ) = 


mL 


E 


ATifck n=1 nsinh(^) 


kdJ 0 (2nkd ) -|— Ji{2nkd) — kdJ 2 (2nkd )) , (B7) 

n 


and the shell energy: 


1 In the Refs. [18,11]] the factor “l/7r” has been missed in the expression for the shell energy and 
level density of the rod-like phase. 
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ti 2 k 2 T 00 


8rmr d ~j n 2 sinh(^) 


-J2(2nkpd). 


(B8) 


d 


In the limit — >> 1, using the asymptotic formula for a Bessel function and expressing the 

R 

Fermi momentum as a function of the density, one gets: 


9c d (k,d ) RJ t ms ( 2kd ~ j) + 0(l/d 312 ), 


(B9) 


and for the shell energy: 


f, 2 T R / Q \ l/ 2 n 1 / 2 7T 

KhtlAd) « (-) ^575 ws( 2 nk F d - -) + 0 (l/d T/2 ). 


(BIO) 


In the case of two slabs the semiclassical approximation is no longer applicable since 
the classical orbits in this case are parabolic which implies that A = 0 and the expression 
for gc diverges. On the other hand, the problem can be solved analytically yielding an 
exact solution for the shell energy ||T8[. In the following, however we derive an expression 
in a different way, which has the advantage that the final formula has a form similar to the 
semiclassical expressions for spheres and rods. 

We consider first the fluctuating level density of a one-dimensional Fermi gas confined 
in a box of length d with Dirichlet boundary conditions. The total level density in this case 
can be calculated exactly and has the form 


(Bll) 


md 00 

g(k,d) = 2 (1 + 2E cos (2 nkd)). 

n kn n=1 

The second term corresponds to the fluctuating part of the level density: 

2 md 00 


9c dim %d) = 


h~kn 


y, cos(2 nkd). 


(B12) 


n =1 


Using an analogous approach as in the case of rods, one can calculate the fluctuating part 
of the particle number as a function of the momentum: 


= l 9c im (K,i)d 

= f-l 2 -£ 

1 27 t) 2d 2 z - 


'rfk f 

2 m 


! U2_|_£.2^ 1.2 _b.2 
K'x' ^z 


dk x dk y 


n=1 


—2nkpdcos(2nkpd) + sin(2 nkpd) 

q : 

n 6 


(B13) 


where k z is the momentum component perpendicular to slabs. Hence the fluctuating level 
density for slabs is given by: 


9 ?‘( M) 


dNpj ab m 
dk h 2 k 


2m 



sin(2 nkd) 
n 


The shell energy calculated using expression (|A3|) gives: 


(B14) 
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TT'slab _ 

^ shell 


n 2 L 2 d ~ 1 

32mn 2 ( nd ) 5 


(6nkfdcos(2nkpd) + 4 {nkpd) 2 sin (2nkpd) — 3sin(2 nk F d)). (B15) 


Note that unlike in the previous cases the truncation of the summation for large separation 
is not justified since all terms have the same dependence on the distance. Nevertheless, the 
gross shell structure is still represented by the n = 1-term and thus for comparison with 
previous expressions we consider the leading term n — 1: 

E$$,(d) *-(-) P — sin(2 k F d) + 0(1//). (B16) 


and for the fluctuating level density: 


1 TY) T ^ 77 * 

d )^~2 cos(2 kd - -) + 0(1// 


(B17) 


The expressions ( |B3| . |B10| . [1 11 (i| ) and ( [1 >4] . [1 Ttlfl 11 7p can be gathered into a single formula 
for the shell energy and the fluctuating part of the level density, respectively: 


Elheii(d) 


tfUR 2 - 1 

8 m 


Q \ ( \ 2+1 

3\ 6 (p) 


7r 


d^r 


,7T 


— cos ( 2k F d — i— ) + 0(l/d 5 2 ), 


(B18) 


where i 


i+2 


Q i (e d) « l UR 1 (—) 4 
9c[,) 4 d 1_i / 2 \TT 2 h 2 


0,1, 2 correspond to the spherical, 


cos ^2kd — + 0(l/d 3 2 ), 

rod and slab case, respectively and e 


(B19) 

n 2 k 2 

2771 
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FIGURES 


FIG. 1. The total energy density (jlj) of the npe matter confined within the cubic box of length 
d as a function of the proton fraction. The SLy4 force has been used for the nucleon-nucleon 
interaction. All densities are specified in fm -3 . 

FIG. 2. The total energy density (jl|) of the npe matter as a function of the proton quadrupole 
moment Q p = Qfo (middle subfigures). The integrated proton and neutron densities (see text for 
definition) corresponding to nuclear configurations indicated by arrows are shown in the lower and 
upper subfigures. 


FIG. 3. The same as in the Fig.2 but for different densities. 


FIG. 4. The energy density difference A E/V between nuclear phases as a function of the total 
density. Solid curve denotes the difference between the spherical and rod-like phase. Dotted curve 
denotes the difference between the spherical and slab-like phase. Smaller subfigures show the 
energy density of the npe matter as a function of the proton quadrupole moment for four different 
densities. Parameter d denotes the length of the cubic box. 


FIG. 5. The gross structure of the interaction energy (shell energy) as a function of a distance 
between two impenetrable obstacles immersed in the neutron gas at the density p = 0.05 fm~ 3 , 
calculated in the semiclassical approximation. The radius of the spherical obstacle (a) R = 7 fm, 
the radius and length of the cylindrical obstacle (b) is equal to R = 4 fm and L = 28.58 fm, 
respectively. The width and length of the planar obstacle (c) is equal to R = 7 fm and L = 14.32 
fm. The sizes were adjusted to keep the volume of the obstacles constant. 
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